Hybrid molecular-continuum fluid dynamics 
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We describe recent developments in the hybrid atomistic/continuum modelling of dense fluids. 
We discuss the general implementation of mass, momentum and energy transfers between a region 
described by molecular dynamics and the neighbouring domain described by the Navier-Stokes 
equations for unsteady flow. 
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O |. The flow of complex fluids near interfaces is governed by a subtle interplay between the fast microscopic dynamics 
within a small localised region of the system close to the interface and the slow dynamics in the bulk fluid region. 
This scenario is encountered in a wide variety of applications ranging from nanotechnology (nanofluidics) and other 
industrial processes (such as wetting, droplet formation, critical fluids near heated surfaces or crystal growth from 
a fluid phase) to biological systems (for example, membranes or biomolecules near interfaces). The dynamics of 
these systems depends on the intimate connection of many different spatio-temporal scales: from the nanoscale to 
the microscale and beyond. Realistic simulations of such systems via standard classical molecular dynamics (MD) 
. — , are prohibitive, while continuum fluid dynamics (CFD) cannot describe the important details within the interfacial 
region. In view of this fact, the field of computer simulation is now faced with the need for new techniques, which 
bridge a wider range of time and length scales with the minimum loss of information. A hybrid particle-continuum 
O ■ approach provides a resolution to this dilemma. A hybrid algorithm retains all the atomistic detail within the 
'— 1 \ relevant localized domain and couples this region to the continuum hydrodynamic description of the remainder of the 
system. Indeed hybrid algorithms for liquids can be expected to provide a powerful tool for the fast growing field 
of nanofluidics in micro electro-mechanical systems (MEMS) and our ongoing contributions have been recognized by 
, the nanoscience community (R. Delgado-Buscalioni and Coveney 2003a) as offering a promising simulation technique 
' with nanotechnological applications. 

Hybrid algorithms for solids (Abraham et al. 1998) and gases (Garcia et al 1999) were the first to be fully developed 
in the literature. As expected in most theoretical descriptions of matter, the hybrid description of the liquid state is 
the most challenging one. The general procedure is to connect the particle domain (P) and the continuum domain (C) 
within an overlapping region comprised of two buffers: C— >P and P— >C (see figure 1). Within the P— >C buffer the 
particle dynamics are coarse-grained to extract the boundary conditions for the C-region. The most complicated part 



of any hybrid scheme is the C^P coupling where the microscopic dynamics need to be reconstructed to adhere to 
. — ' the prescriptions given by the continuum variables. Moreover, in doing so the unphysical artifacts thereby introduced 
, should be minimized (following Occam's razor). 

In this paper we provide an overview of the state-of-the-art of the hybrid modelling of liquids. In ij^we start by 
presenting an overview of the hybrid scheme and some preliminary topics such as the inherent constraints on the 
continuum time step and the spatial-grid size. Section fllll discusses several implementations of the temporal-coupling. 
The C— >P coupling scheme is then explained in flIVI for the general case of mass, momentum and energy. We illustrate 
this important part of the scheme by reproducing the three hydrodynamic modes (shear, sound and heat) governing 
the relaxing flows in an infinite medium. Section[V]is devoted to the P— >C coupling, based on a finite volume method 
solving the flow within the C domain. Some comments on the effect of noise on the accuracy of the scheme are made. 
The full method is used in HVII to solve the problem of shear flow driven by oscillatory wall motion in a nano-slot. 
Finally, conclusions and future directions for this research are described in Will 
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II. OVERVIEW 

The domain decomposition deployed in our hybrid scheme is depicted in figure 1. Within domain P the fluid is 
described at the atomistic level via Newtonian dynamics. The position of the N(t) atoms at time t inside P is updated 
each Atp time interval using a standard MD scheme. The present calculations were done with a Lennard-Jones (LJ) 
fluid. Throughout the ongoing discussion all quantities are given in reduced Lennard Jones units: length er, mass 
m, energy e, time (m<T 2 /e) 1 / 2 and temperature e/fcp. We refer to Hoheisel f996, for the estimated physical values 
of the LJ parameters for several substances (as an example, for a simple molecular fluid as N2, a ~ 0.35nm and 
e/k B ^ 100K). 

The rest of the computational domain (C) is described by the Navier-Stokes equations. The fluid variables at C are 
the densities of the conserved quantities for which the equations of motion in conservative form are d^/dt = — V • J$ 
with <I> = {p, pu, pe} and J$ = {pu, pu + II, pue + II • u + q} standing for the mass, momentum and energy fluxes 
respectively. Here p is the density, u the local velocity, e the specific energy, II = PI + t the stress tensor which 
contains the pressure P and the viscous tensor r (for a Newtonian fluid) and q = — kV-T, the heat flux by conduction 
expressed via Fourier's law. These continuum equations may be solved via standard CFD methods. Alternatively, for 
low- Reynolds number flows (Re < O(10)) the equations can be solved analytically (Delgado-Buscalioni & Coveney 
2003b), as is done in the tests presented in illVI 

The kind of information to be transferred in the overlapping region has been the subject of some discussion. The 
first attempts in the literature (Delgado-Buscalioni & Coveney 2003b and references therein) considered the transfer 
of momentum in steady shear flows and proposed a matching procedure based on the continuity of velocity across 
the overlapping region. This sort of coupling strategy may be referred to as "coupling-through-state" . An alternative 
formulation of the information exchange for liquids based on matching the fluxes of conserved quantities (to/from P 
and C) was proposed by Flekk0y et al. 2000. These authors considered steady shear flows with mass transfer. In 
subsequent work by Delgado-Buscalioni & Coveney (2003b) the flux-coupling scheme was generalized to enable transfer 
of mass, energy and momentum (along both transversal and longitudinal directions). Delgado-Buscalioni & Coveney 
2003b also present a comparative study of the coupling-through-fluxes and coupling-through-state schemes for flows 
involving energy transfer (longitudinal waves). It was shown that the coupling-through- fluxes scheme provides the 
correct physical behaviour, while the coupling-through-state scheme does not guarantee positive entropy production. 
Consequently the coupling of fluxes is of central importance in our hybrid scheme (see fllVI and flVjl . 

III. TEMPORAL COUPLING 

In general there are three times involved in the coupling scheme: the MD time-step Atp, the time-step for the 
C-solver Atc(>> Atp) and the averaging time At av , which are presented below as outline two possible strategies for 
merging the time-evolution of C and P. The information transfer (from C— >P and P^C) is updated over each time 
interval, Ate- As stated above, the P— >C coupling consists firstly of a coarse-graining procedure. In particular, for 
any particulate quantity, <E>i, the spatial average over each P— >C cell of volume Vpc {= AAXpc, in figure 1) is defined 
as $(R, t) = T<i£Vpc®i/Npc, where R is the position of the cell in the coarse-grained coordinates and Npc is the 
number of particles inside Vpc- The time average also needs to be local with respect to the coarse-grained dynamics. 
To that end, the microscopic quantities are sampled over a time interval At av which is treated as an independent 
parameter of the simulation: 

(*) (R,t c ) = -L- f tC $(R,t)dt. (1) 

Ate Jt c -At a v 

The magnitudes Ate an d At av are constrained by several physical and numerical prerequisites quoted in Table 1. 

There are essentially two ways to deal with the coupling of time within the hybrid scheme: sequential coupling 
or synchronized coupling. The diagrams in fig. 2 illustrate two possible choices for these time-coupling strategies 
starting from given initial conditions. In the sequential coupling scheme, both P and C are first moved to t = Ate 
using the initial conditions. The C— >P coupling is performed first at t — Ate and the P system is advanced to 
t = 2 Ate ~ 300 Atp. The averaged P-information collected over time interval At av = 2 Ate within the P— >C 
cell is then transferred to the C domain, giving the required boundary condition to advance C towards the same 
time t — 2 Ate- This procedure is suited for serial processing. More refined versions of sequential coupling can be 
constructed to perform averaging over times At av greater than Ate ■ 

In the synchronized coupling scheme both domains advance in time independently until a certain instant at which 
both C— >P and P— >C information transfers are exchanged. This scheme is suitable for parallel processing because the 
P and C domains are being solved concurrently. We note that in this case the averaged information from P transferred 
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Physical condition 


Algcbraiconstraints Eq. 


Local Equilibrium 
Flow resolution 
Accuracy 
Courant condition 


At c > r co i = 0.14p _1 T- 1/2 Ax > A = 0.2P- 1 C.l 
At av < O(O.l) T flow I*" 1 d$/dx|Aa; < 1 C.2 
VpcAt„„ > r/( 7 2 »7) C.3 
At c < Ax/{2u flow ) C.4 



TABLE I: Constrains on the coarse grained time and length scales within our hybrid MD-CFD scheme. Condition C.l ensures 
the local thermodynamical equilibrium at the averaging region: the coarse-graining time Atav and grid-spacing Arc needs to be 
larger than the collision time r co ; and the mean free path A, respectively. In C.l r co ; and A are estimated by the hard-sphere 
approximation. Condition C.2 is needed to resolve the fastest flow characteristic time Tfi ow and the spatial variation of any 
physical variable $ over the control cell |$ _ d<& / dx\Ax. Depending on the flow behaviour, in C.2, Tfi ow may stand for the 
period of the oscillatory flow / _1 or for the diffusive time L 2 x /v, etc. The accuracy condition in eq.(C.3) ensures that the signal- 
to-noise ratio of the transversal momentum flux in a flow with shear rate 7 is greater than one (similar kind of relationships 
can be derived for the longitudinal momentum and energy fluxes). The conditions C.l, C.2 and C.3 are applied within the 
P^C cell. The last condition, C.4, ensures the stability of the numerical (explicit) scheme used for time- integration of the C 
flow. The characteristic velocity of the flow (on one grid space) is denoted by Ufi ow - 

at any of these times is obtained during the previous time interval At av . This fact introduces a delay of 0(A£ o „/2) in 
the C flow. Hence, it is important to ensure that At av is about O(10 _1 ) times smaller than the fastest physical time 
of the flow process (see Table 1). 

IV. CONTINUUM-TO-PARTICLE COUPLING AND ITS VALIDATION 

The generalized forces arising from fluxes of mass, momentum and energy measured from the C flow are to be 
injected into the particle system at the C^P cell. Table 2 summarizes how each flux contribution arising within C is 
translated into the P domain. 

Mass continuity is ensured by inserting or extracting particles at a rate given by eq. (T.l) in Table 2. The 
convection of momentum is determined by the product of the rate of particle insertion s and the average velocity of 
the incoming/outgoing particles (v'). By injecting eq. (T.l) into eq. (T.2) it is easily seen that convection balance 
requires (v') = u. New particles arc therefore introduced with velocities sampled from a Maxwellian distribution 
at temperature T and mean velocity u. On the other hand, the local equilibrium (v) = u ensures that the average 
velocity of any extracted particles is equal to that of the continuum prescription 

Viscous and pressure forces are introduced via external forces acting on the particles at P^C. An important issue 
is to decide how to distribute the overall force in eq. (T.3), J2 Npc -p{ ext \ over the individual particles. We refer to 
Flekk0y et al. (2000) and Delgado-Buscalioni & Coveney (2003b) for a full discussion. Although in general the force 
to be felt by each particle i within the P^C cell can be distributed according to the particle positions (see Flekk0y 
et al. 2000), we have adopted a flat distribution Ff xt = All ■ n/Npc because it provides by construction, a correct 
rate of energy dissipation in eq. (T.5) (sec Delgado-Buscalioni & Coveney 2003b). Using eq. (T.l) it is seen that 
the balance of advected energy in eq. (T.4) implies (e') = e. The energy of each particle is composed of kinetic and 
potential parts, et — vf/2 + "4>i({r} ). The specific energy of the continuum is e = u 2 /2 + 3kT/(2m) + 4> (here (j) is 
the excess potential energy). The balance of kinetic energy {{v') 2 ) — u 2 /2 + 3fcT/(2ra) is ensured by inserting the 
new particles with the proper Maxwellian distribution. The balance of the potential energy requires a more difficult 

condition (i>{{r} N )^ — 4> to be satisfied. When inserting a new particle, this involves finding a precise location within 

the C^P cell with the desired potential energy. To solve this problem in a fast and effective way we have constructed 
an algorithm for particle insertion called usher (Delgado-Buscalioni & Coveney 2003c). In order to find the site 
with the desired energy within the complex potential energy landscape, the usher algorithm uses a variation of the 
steepest descent algorithm including an adaptable displacement. For densities within the range 1 p = [0.4 — 0.8], the 
usher scheme needs around 8 — 30 iterations, each one involving the evaluation of a single-force. The usher algorithm 
can be also applied in other problems involving particle insertion, such as grand-canonical molecular dynamics. 

Finally, eq. (T.6) in Table 2 determines the rate of heat transfer into P by conduction. This energy can be injected by 
reproducing a non-isothermal environment within the C— >P cell. To that end we have implemented a set of (typically 
2-3) Nose-Hoover thermostats (NHT) separated by a distance d with temperatures differing by AT = [VT-n] d, where 
VT is the C-temperature gradient at C^P. 

The decay of transversal and longitudinal waves is an excellent test for the validity of our proposed the C^P 
coupling as they comprise the whole set of hydrodynamic modes: shear, sound and heat waves. For these tests we 
implemented a set-up consisting of a P region of length L x (with periodic boundary conditions in y and z directions) 
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Conserved 
quantity 


Fluxes 


P 


«- c 


cq. 


Mass 




ms 


= Apu ■ n 


(T.l) 


Momentum 


Convection 


ms (v ) 


= Apuu ■ n 


(T.2) 




Stress 


\L, F i / 


A 1 1 , . 

= All ■ n 


(T.3) 


Energy 


Advection 


ms (e') 


= Apue ■ n 


(T.4) 




Dissipation 




= ATI u n 


(T.5) 




Conduction 




= Aq ■ n 


(T.6) 



TABLE II: The balance of mass, momentum and energy fluxes at each C — > P cell. The fluxes measured within C (third 
column) are imposed into P via the expressions given at the second column. The cell's surface is A, and the surface vector n 
points outwards (fig. I). The mass rate is s(t) (s>0 for inserted and s<0 for removed particles). The velocity and energy of 
the inserted/removed particles are v' and e' respectively. The external force and heat flux inserted within the C^P cell are 
J2 Npc F? xt and <J^*>, respectively. 

surrounded by two C domains. We initially imposed on the P system a sinusoidal (x- or y-) velocity profile along 
the x direction. By extracting the initial amplitudes of the spatial Fourier components of all the hydrodynamic 
quantities it is then possible to trace the entire time-evolution of the relaxing flow using linear hydrodynamics. In 
particular, this permits us to calculate at any time the generalized forces to be inserted within the C— >P cell. The 
time evolution of the spatial Fourier components of the P-variables is finally compared with the analytical expressions. 
Such kinds of comparisons are shown in figures 3 and 4, for the case of a relaxing shear wave and a longitudinal wave, 
respectively. The excellent agreement obtained indicates that the C— >P coupling protocol can be used for capturing 
fast and low-amplitude flows, such as those governed by sound, shear or heat waves. 

The entropy perturbation, shown in fig. 5, was calculated from the temperature and density perturbative field. The 
results clearly show that using only one thermostat per C— >P cell (denoted by l-NHTcp, in figs. 3 and 4) leads to 
negative entropy production. The pure exponential decay of heat due to diffusion is only recovered when the correct 
(averaged) heat flux is connected to each C— *P cell; in figs. 4 and 5 we present a result with two thermostats per cell 
(2-NHTcp). This result confirms that the coupling-through-fluxes scheme is the correct matching procedure. 

V. PARTICLE-TO-CONTINUUM COUPLING: FINITE VOLUMES AND FLUCTUATIONS 

Within the P— >C cells the information coming from the particle dynamics is coarse-grained to provide boundary 
conditions at the "upper" C-level. In i[n]we introduced the averages needed to produce such information. At the 
P^C interface the C region receives the averaged particle-fluxes as open-flux (von Neumann) boundary conditions. 
The averaged mass, momentum and energy particle-fluxes through the P^C interface are constructed as follows, 

pu-n PC = (s^mvA • npc (2) 

V pc \ / 

n • np C = ^ femViVi - ±E?f '•<•/■•'<;)) ' n ?c (3) 

q • = \ pc ^ (s£°e iVi - ^, V r ••- ; v,F, ; ) ^ • n PC (4) 

where Npc is the number of particles inside the P^C cell of volume Vpc and npc is the surface vector shown in fig. 
1. 

A. Hybrid finite volume: boundary conditions 

Let us now illustrate how these fluxes can be injected into the C domain in the framework of the finite volumes 
method (Patankar 1980). The finite volumes method is ideally suited to our scheme because it exactly balances the 
fluxes across the computational cells. Its principle is simple. Briefly, the computational domain (C) is divided into 
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cells of volume Vi whose size and location is given by the nodes of a specified mesh, {R;}, I = {1, M c }. Integrating 
the conservation equation d$/dt = — V • J$ over each computational cell (say the cell H in fig. 1) one obtains, 

/ 

where Af stands for the area of the face / and is the outwards normal surface vector. The volume integral of the 
transient term of the conservation equation has been approximated by Vjj times the explicit time derivative of the 
value of the integrand at the cell centre, halfway between the surfaces: ph^h- Equation (JjjJ yields a set of ordinary 
differential equations (ODE's) involving the flow variables at each cell face, /. The set of equations is closed for the 
flow variables at the cell centre by expressing the fluxes at the interfaces J$ ,/ in terms of differences of flow variables 
at neighbouring cell centres, via the constitutive relations. 

Let us consider the momentum flux balance for the low Reynolds number flow of an incompressible and isothermal 
fluid driven by diffusion of y-velocity along x direction: u = u(x)j. In this case J • n = Pi — r)(du/dx)j, where 
the surface vector of the P-^C surface is n — i. Let us consider an isobaric environment and restrict ourselves to 
the transfer of transversal (y) momentum, governed by the momentum flux J = J • j = —rjj and the shear rate 
7 = du/dx. Integrating along the cell H (see fig. 1), using a first order space discretisation of the stress (e.g. 
J w = —r)(uH ~ uw) I 'Ax) and an explicit time integration scheme, one obtains 

u H = u H (l - 2r) + ru E + ruw, (6) 

where the subscripts H denote the set of cell centres H = {1,M}, and the symbols E (east) and W (west) denote 
variables measured: x = E{= H + 1) and x = W(= H — 1). The time instant is denoted by uh — u(xH,t) and 
U H = u(xh, t + Ate) and r = vAt/(Ax 2 ). with v — r//p the kinematic viscosity. In order to guarantee the numerical 
stability of the explicit scheme in eq. ©, the size of the (smallest) control cell inside the C region Ax and the time 
step At are related through r < 1/2, which corresponds to the grid-diffusive- velocity Ufi ow — v/ Ax in the Courant 
condition C.4 of Table 1. In solving eq. JfjJ) we used a uniform grid with a typical value of Ax ~ 0.5. 

In order to impose the boundary condition one needs to determine the velocity within the outer cells: at the 
rightmost x = xm+i = L x and at the leftmost boundary (inside the P— >C cell, see fig. 1) x — xq = lc — Ax/2. At 
%M+i = L x there is a rigid wall which moves at a velocity u wa u(t) and provides the Dirichlet boundary condition 
UM+i — u waii(t)- The hybrid formulation is applied at the left boundary xq = lc — Ax/2. To evaluate the outer 
velocity uw = uo we impose the balance of momentum flux across the w surface at x = lc- This means that the 
continuum flux evaluated at x = w is made equal to the corresponding averaged particle flux (j) w = —t](uh~uw)/Ax. 
The outer velocity to be inserted in eq. I© is then uw — uh + {j)wAx/rj. The velocity uh is evaluated as a linear 
combination of the continuum uh{= u\) and the average particle velocity (v)h at x\ = lc + Ax/2: 

uh = (1 - o)uh + o>(v)h- (7) 
By inserting eq. J7J into eq. © one obtains the velocity at the boundary cell 

Uh = u H (1 - r) + r u E + i-^i + ar ((v) H - u H ) ■ (8) 

pAxH 

The reason for the choice of uh in eq. now becomes clear. It introduces the last term on the right hand side 
of eq. which acts as a forcing term ensuring velocity continuity by gently driving the continuum velocity to the 
corresponding particle average ur — {v)h- The strength of the velocity coupling is maximal when a — 1 and is absent 
if a = 0. The idea of using a hybrid gradient (arising for any a ^ in eq. J7J) arose from the outcome of calculations 
performed at very low shear rates (7 < 10~ 2 ). Using a — one obtains a velocity discontinuity at P^C which is 
of the same order of magnitude as the fluctuations of the mean instantaneous velocity within the overlapping region. 
At low shear rates this means substantial relative differences in the C and P velocities, < v >h —uh)/uh ~ O(l). 
This problem is solved by introducing a small velocity coupling in the continuum scheme, with a small value of 
a G [0.2,0.5], which drives the continuum velocity to the average particle velocity in a time of 0[Ax 2 / \va)]. To check 
any influence of the velocity coupling term in eq. JSJ on the flux balance, we performed simulations of the Couctte 
flow at different shear rates and compared its average over time with the time averaged momentum particle flux. The 
results showed that, in average, the velocity coupling term is vanishingly small so it does not introduce any extra flux 
in the coarse-grained time scale. 



B. The effect of fluctuations: shear stress 



In our scheme, the fluctuating nature of the fluxes introduced into the C region at P^C imposes a limitation on 
our ability to resolve the flow field, as also arises in experiments and full MD simulations. This limit is determined by 
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signal-to-noise ratio becoming smaller than one. A theoretical expression for the amplitude of the stress fluctuations 
can be obtained (Delgado-Buscalioni et al. 2003), providing a relationship between the signal-to-noise ratio and the 
coarse-grained time and space scales At av and Vpc- Table 1 contains the condition to ensure an averaged shear force 
larger than its variance. It is clear that in weak steady flows it is always possible to increase the signal-to-noise ratio 
by enlarging At av . Nevertheless, in a general space and time-dependent flow, the sizes of the averaging windows in 
space and time (Vpc and At av ) are bounded above by the minimum wavelength and characteristic time which need 
to be treated within the flow. Such requirements on spatial and temporal flow resolution are also quoted in Table 1. 



In order to test the applicability of the full hybrid scheme under unsteady flows, we have considered the flow of an 
incompressible and isothermal fluid between two parallel walls in relative oscillatory motion. This set-up is widely 
used to investigate the reological properties of complex fluids attached to surfaces, as polymer brushes (see CM. 
Wijmans & B. Smit (2002) for a recent review). These systems are good examples of the sort of applications of the 
hybrid scheme, which can treat the complex fluid region by MD and the momentum transfer from the bulk by CFD. 
A similar set-up can be also used in the simualation of nanotechnological process. For instance, Stroock et al. (2002) 
showed that the mixing of solutions in low Reynolds number flows in microchannels can be enhanced by introducing 
bas-relief nano-structures on the floor of the slot. In our test flow, the simulation domain is < x < L x and it is 
periodic along y and z directions. The particle domain occupies the region x < lp, and it includes the LJ liquid 
and the atomistic wall composed of two layers LJ particles at x < 0. The continuum domain comprises the region 
x G [lc, L x ]. The sizes of the simulation domains were within the nanoscale L x ~ 50<r, and lp ~ 15<7, while the width 
of the overlapping region, lp — lc, was set to arround 5cr. The flow is uniquely driven by the oscillatory motion of the 
x = L x wall along the y direction, meaning that the mean pressure is constant throughout the domain and there are 
no transfers of mean energy or mass in the x direction (perpendicular to the P— >C surface). Therefore the mean flow 
carries transversal momentum by diffusion only, and the equation of motion for the y-velocity is du/dt = vd 2 u/ 'dx 1 , 
with boundary conditions u(0,t) = and u{L,t) = u wa u{t) = it max sbx(u>t). This equation can be solved analytically 
(H. Schliting 1958; CM. Wijmans & B. Smit 2002). The flow profile has a maximum amplitude at the moving wall 
and the momentum introduced by its motion penetrates into a fluid layer of width <5 ~ tJiw/ f. Beyond this layer 
the flow amplitude tends to zero diffusively as it approaches the other wall held at rest. Therefore, the maximum 
shear rate attained inside the momentum layer is of order 7 ~ u max /S. Inserting this relation into the signal-to-noise 
condition (C.3 in Table 1), we find 



Equation means that in order to attain a signal-to-noise ratio larger than one, the mean kinetic energy per unit 
volume of the flow integrated over the averaging time At av needs to be larger than the corresponding energy due to 
fluctuations over the period of the mean flow. It is important to mention that at low enough frequencies (/ > v/ L 2 ), 
there is sufficient time for momentum to be spread by diffusion over the whole domain. In such situations the correct 
condition is given by the signal-to-noise condition (C.3 in Table 1) with 7 ~ u max /L 2 . 

As indicated by condition C.2 in Table 1, in order to solve for the temporal variation of the flow it is required that 



At av f > O(0.1). Inserting this condition into eq. © one obtains u max > 5 ( k ^ T ) . For k B T = 1.0, p = 0.8 



and Vpc = O(100) the above inequality yields u ma x > 0.5. We performed oscillatory shear simulations for values 
of u max above, close to and below the threshold given by eq. ©. As shown in fig. 6a, calculations made at large 
flow amplitudes are in excellent agreement with the analytical solution. In figure 6b we present results for the same 
density and temperature (p — 0.8 and T = I) and a wall velocity u max = 0.5 right at the accuracy limit predicted by 
@. The averaging time was chosen to be At av — 10. As shown by the instantaneous velocity within the P^C cell, 
the noise amplitude is nearly equal to the flow amplitude and its time-averaged value shows traces of fluctuations. 
Figure 6c corresponds to the same velocity and density but at a larger temperature T = 4. This case is below the 
accuracy limit (given by C.3 in Table 1) where forces arising from thermal fluctuations dominate the hydrodynamic 



VI. OSCILLATORY WALL FLOW 




(9) 




ones. 



VII. CONCLUSIONS AND FUTURE DIRECTIONS 



We have presented a hybrid continuum-particle scheme for moderate-to-large fluid densities which takes into account 
mass, momentum and energy exchange between a domain described by discrete particle Newtonian molecular dynamics 



7 



(P) and an interfacing domain described by continuum fluid dynamics (C). The coupling scheme is applied within 
an overlapping region comprised of two sub-cells where the two-way exchange of information is performed: C^P 
and P— >C. We have shown that the coupling-through-variables scheme (which simply ensures continuity of variables 
within the overlapping region) is not sufficient to guarantee positive entropy production. However, by generalizing 
the coupling-through- fluxes scheme proposed by Flekk0y et al., 2000 to energy and mass transfer we find that the 
correct decay of shear, sound and heat waves is obtained. 

We are now deploying the present scheme to study the dynamics of a tethered polymer under shear flow. The 
polymer and its local environment are treated via MD, while the shear flow imposed on the outer domain is treated 
via the finite volume CFD method. In the future, we plan to apply our hybrid scheme to the study of membrane 
dynamics. 

Enhancements to the present hybrid algorithm are under investigation. In the scheme described here the energy 
flux balance is ensured only over time averages. We are currently studying alternative schemes which exactly balance 
this flux. From a numerical standpoint, we plan to implement the P— >C coupling in conjunction with a finite volume 
CFD solver in 3D. 

Also, the present scheme can be easily adapted to couple molecular dynamics with another mesoscopic scheme that 
takes into account hydrodynamic fluctuations. This sort of hybrid scheme could be used in applications where the 
fluctuations are relevant (microfluidics, fluids near critical point, etc.). An important condition for the interfacing 
mesoscopic scheme is that it needs to be fully consistent with thermodynamics. Also important is that the transport 
coefficients of the mesoscopic model should be adjustable to represent the correct coarse-grained dynamics of the 
selected working fluid. Natural candidates are the Lagrangian schemes involving Voronoi tesselation (Flekk0y et 
al. 2000a) or the Smooth Particle Dynamics model and related mesoscopic techniques (Espahol 2003). The lattice 
Boltzmann (LB) method is another possible candidate to interface with the MD domain. This model has been already 
used in multiscale modelling (Succi et al 2001). Nevertheless, the problem with LB methods at present is that there is 
no truly reliable thermohydrodynamic model other than for single phase flow. Energy conservation remains unsolved 
and most models are athermal; even the thermohydrodynamic lattice-BGK models for the ideal gas are vastly over- 
determined and get the temperature dependence of the viscosity wrong (Boghosian and Coveney 1998). Therefore 
the hybrid scheme proposed here could only be interfaced with the lattice Boltzmann model in certain applications 
involving isothermal and incompressible single phase flows. 

A longer term goal of this research is to develop a flexible, componentized, hybrid coupling environment into which 
any molecular dynamics and any continuum fluid dynamics codes may be inserted. This will require considera- 
tion of electrostatic forces and, therefore, an additional conserved quantity, the electric charge, whose flux coupling 
will requires use of Poisson-Boltzmann solvers. Moreover, such multiscale hybrid schemes are attractive candidates 
for efficient deployment on computational grids, a feature now under investigation with the RealityGrid project 
(www.realitygrid.org). 
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P->C 

FIG. 1: The domain decomposition of the hybrid scheme: (a) displays P and C regions separatedly. The shaded region 
represents the overlapping domain comprised by a 2D array of C— >P and P^C cells where the exchange of microscopic and 
macroscopic information is carried out. The surface area of each cell is A. (b) shows the P^C region in more detail and the 
neighbouring control cells pertaining to the finite volume discretization of the C region. In this one-dimensional example, the 
width of the P^C cell is Axpc and its volume is Vpc = AAxpc- 
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FIG. 2: Two possible time coupling strategies in a particle-continuum hybrid scheme: (a) synchronized coupling and (b) 
sequential coupling. Bold arrows indicate the direction of the information transfer. The time average of the P variables is 
performed during the time interval At a v by n 3 samplings separated in time by St s . Ate and Atp are the continuum time step 
and the MD time step respectively. 
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FIG. 3: The Fourier components of the transversal and longitudinal velocity perturbation. In the notation Vy ; s indicates 
the sinusoidal component and c cosinusoidal and n the wavenumber k n — nko. The transversal wave has ko = 0.35, the size 
of the P region was and L x = 20 (in LJ-units) and the temperature was T = 2.5; while for the longitudinal wave fco = 0.168, 
L x = 40 and T = 3.5. In both cases p — 0.53 The autocorrelation of the velocity is also shown. In all graphs the dashed lines 
are the analytic solution from linear hydrodynamics. Reproduced from R. Delgado-Buscalioni & Coveney 2003b with permision 
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FIG. 4: The dominant Fourier mode of the various thermodynamic variables in the decay of the same longitudinal wave 
shown in figure 3. Comparison is made between a calculation with tow Nose-Hoover thermostats per C^ P cell (2-NHTcp) 
and another using only one thermostat (1-NHTcp). Dashed lines are the analytical hydrodynamic solution. The entropy 
production from these two simulations is shown in figure 5, only the one with two thermostats yields th correct physical 
behaviour. Reproduced from R. Delgado-Buscalioni & Coveney 2003b with permision 
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FIG. 5: The main Fourier mode of the entropy density (as a product with the mean temperature) — (^Q^' 1 ^ time-averaged 

along Atav = 1.0. The result comes from the same longitudinal wave shown in figs. 3 and 4. Comparison is made between 
a flux-coupling scheme (using 2-NHTcp) and the coupling-state scheme using 1-NHTcp (cf. fig. 4). The latter violates the 
second law of thermodynamics. The dashed line is the analytical hydrodynamic result. Reproduced from R. Delgado-Buscalioni 
& Coveney 2003b with permision. 
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FIG. 6: Mean molecular velocities within the overlapping region for several oscillatory-wall shear flows applied to a LJ fluid; 
^max is the maximum wall velocity and f its frequency, (a) Flow corresponds to tt m ax = 10, / = 0.01 and T = 1.0; we plot the 
instantaneous particle velocity at P^C and the time-averaged particle velocity (along Atav = 1) at C-^P; (b) corresponds to 
M ml x = 0.5, / = 0.01 and T = 1.0; (c) to u max = 0.5, / = 0.01 and T = 4.0. In all cases p = 0.8, the extent of the periodic 
directions are Ly — L z — 9, while Vpc — Ax L y L z — 178. In (b) and (c) we show the P — >C mean velocity (instantaneous and 
time-averaged velocity with At av — 10); dashed lines are the analytical hydrodynamic solutions of the imposed shear flows. All 
quantities are given in reduced LJ units. 



